! Interpolation of mesh function to the Lagrangean particle mesh
! Case of integer nodes (cell borders) 
! F_{i,k} = F(x_i, Y_k) 
      function inpol1(x,y,hx,hy,f, im, km) 
      real f(im,km),x,y,hx,hy,s1,s2 
      integer im,km,i,k 
      s1=x/hx 
      i=int(s1+ 1) 
      s1=s1-i+1 
      s2=y/hy 
      k=int(s2+1)
      s2=s2-k+1
      inpol1=(1.-s1)*((1.-s2)*f(i,k)+s2*f(i,k+ 1))
     &       + s1*((1.-s2)*f(i+ 1,k)+s2*f(i+1,k+1))
      return 
      end 

! Case of half-integer nodes (cell centers)
! F_{i-1/2, k-1/2} =  F(hx*(i- 1/2), hy*(k- 1/2))
      function inpol2(x,y,hx,hy,f, im, km) 
      real f(im,km),x,y,hx,hy,s1,s2
      integer im,km,i,k
      s1=x/hx
      i=int(s1+1.5)
      s1=s1-i+1.5
      s2=y/hy
      k=int(s2+1.5)
      s2=s2-k+1.5
      inpol2=(1.-s1)*((1.-s2)*f(i,k)+s2*f(i,k+ 1))
     &          + s1*((1.-s2)*f(i+1,k)+s2*f(i+ 1,k+ 1))
      return 
      end 

! Interpolation of a scalar mesh  function from Lagrangean grid
! to the nodes of Eulerian  mesh

! Case of integer nodes (cell borders) 
! rho_{i,k} = rho(x_i, Y_k) 
      subroutine density1 (x,y,mas,jm,ro,im, km,hx,hy,ng1,ng2)
      integer jm,im,km,ng1,ng2,j,i,k
      real x(jm),y(jm),mas(jm)
      real ro(im+ 1,km+ 1)
      real hx,hy,s1,s2
      do 1 i=1,im+1 
      do 1 k=1,km+1
      ro(i,k)=0.
    1 continue 
      do 2 j=1,jm 
      s1=x(j)/hx 
      i=int(s1+ 1) 
      s1=s1-i+1 
      s2=y(j)/hy 
      k=int(s2+ 1) 
      s2=s2-k+1 
      ro(i,k)=ro(i,k)+mas(j)*(1.-s1)*(1.-s2) 
      ro(i,k+ 1)=ro(i,k+ 1) +mas(j)*(1.-s1)*s2 
      ro(i+ 1,k)=ro(i+ 1,k)+mas(j)*s1*(1.-s2) 
      ro(i+1,k+1)=ro(i+1,k+1)+mas(j)*s1*s2
    2 continue
      if(ng1.eq. 1) then 
      do 3 k=1,km+1
      ro(1,k)=2.*ro(1,k)
      ro(im+ 1,k)= 2. *ro(im+ 1 ,k)
    3 continue 
      endif
      if(ng1.eq.2) then 
      do 4 k=1,km+1 
      ro(1,k)=ro(1,k)+ro(im+1,k) 
      ro(im+1,k)=ro(1,k) 
    4 continue 
      endif 
      if(ng2.eq. 1) then 
      do 5 i=1,im+1 
      ro(i,1)=2.*ro(i,1) 
      ro(i,km+ 1)=2.*ro(i,km+ 1) 
    5 continue 
      endif 
      if(ng2.eq.2) then 
      do 6 i=1,im+1 
      ro(i, 1)=ro(i, 1) +ro(i,km+ 1) 
      ro(i,km+ 1)=ro(i, 1) 
    6 continue 
      endif 
      do 7 i=1,im+1 
      do 7 k=1,km+1 
      ro(i,k)=ro(i,k)/(hx*hy) 
    7 continue 
      return 
      end 


! Case of half-integer nodes (cell centers)
! rho_{i-1/2, k-1/2} =  rho(hx*(i- 1/2), hy*(k- 1/2))
      subroutine density2(x,y,mas,jm,ro,im,km,hx,hy,ng1,ng2) 
      integer jm, im,km, ng 1,ng2,j,i,k 
      real x(jm),y(jm),mas(jm) 
      real ro(im+2,km+2) 
      real hx,hy,s1,s2 
      do 1 i=1,im+2 
      do 1 k=1,km+2 
      ro(i,k) =0. 
    1 continue 
      do 2 j=1,jm 
      s1=x(j)/hx 
      i=iat(s1 + 1.5) 
      s1=s1-i+1.5 
      s2=y(j)/hy 
      k=int(s2+1.5) 
      s2=s2-k+1.5 
      ro(i,k)=ro(i,k) +mas(j)*(1.-s 1)*(1.-s2) 
      ro(i,k+1)=ro(i,k+1)+mas(j)*(1.-s1)*s2 
      ro(i+1,k)=ro(i+1,k)+mas(j)*s1*(1.-s2) 
      ro(i+1,k+1)=ro(i+1,k+1)+mas(j)*s1*s2 
    2 continue 
      if(ng1.eq. 1) then 
      do 3 k=1,km+2 
      ro(2,k)=ro(1,k)+ro(2,k) 
      ro(1,k)=ro(2,k) 
      ro(im+ 1,k)=ro(im+ 1,k)+ro(im+2,k) 
      ro(im+2,k)=ro(im+ 1,k) 
    3 continue 
      endif 
      if(ng1.eq.2) then 
      do 4 k=1,km+2 
      ro(2,k)=ro(2,k)+ro(im+2,k) 
      ro(im+2,k)=ro(2,k) 
      ro(im+ 1,k)=ro(im+ 1,k)+ro(1,k) 
      ro(1,k)=ro(im+1,k) 
    4 continue 
      endif 
      if(ng2.eq. 1) then 
      do 5 i=1,im+2 
      ro(i,2)=ro(i, 1)+ro(i,2) 
      ro(i, 1)=to(i,2) 
      ro(i,km+ 1) =ro(i,km+ 1)+ro(i,km+2) 
      ro(i,km+2)=ro(i,km+ 1) 
    5 continue 
      endif 
      if(ng2.eq.2) then 
      do 6 i=1,im+2 
      ro(i,2)=ro(i,2)+ro(i,km+2) 
      ro(i,km+2)=ro(i,2) 
      ro(i,km+ 1) = ro(i,km+ 1)+ro(i, 1) 
      ro(i, 1)=ro(i,km+ 1) 
    6 continue 
      endif 
      do 7 i=1,im+2 
      do 7 k=1,km+2 
      ro(i,k)=ro(i,k)/(hx*hy) 
    7 continue 
      return 
      end 